Experimentos en Redes
Valentina Cardona & Natalia Perdomo
Email vcardonas@unal.edu.co & naperdomol@unal.edu.co
GitHub https://github.com/vcardonas
Rpubs https://rpubs.com/vcardonas
1 Introducción
En todas las ciencias existe un interés generalizado en evaluar el efecto de tratamientos o intervenciones de diversos tipos. Podría decirse que el ejemplo prototípico de este tipo de estudios es el ensayo aleatorio (en inglés, randomized controlled trial, RTC), también llamado como \[Prueba\ A/B\] Este diseño experimental convencional se caracteriza por dos aspectos principales:
La noción de dos grupos a comparar (tratamiento y control, o A y B).
La asignación aleatoria de individuos a dichos grupos.
Sin embargo, en estos experimentos suele asumirse que el tratamiento dado a un individuo no influye en los resultados de otros individuos (no interferencia).
Naturalmente, al realizar experimentos en redes no se puede descartar de manera realista la interferencia.
2 Modelo Neyman-Rubin: Resultados Potenciales
Supongamos que queremos evaluar la efectividad de un tratamiento en un conjunto finito compuesto por \(N\) unidades experimentales. La dificultad de esta tarea radica en que no podemos medir múltiples resultados simultáneamente un individuo dado, aunque deseamos evaluar la diferencia en su variable de resultado bajo diferentes opciones de tratamiento.
El modelo de inferencia causal aborda este problema al incorporar en la medición del impacto la idea de todos los resultados posibles bajo el diseño experimental, es decir, incluyendo tanto los resultados observados como los contrafactuales (potenciales).
| Individuo | \(Z_i = 0\) | \(Z_i = 1\) | \(Z_i = 0\) | \(Z_i = 1\) |
|---|---|---|---|---|
| 1 | \(\mathscr{O}_1(0)\) | \(\mathscr{O}_1(1)\) | \(\mathscr{O}_1(0)\) | \(?\) |
| 2 | \(\mathscr{O}_2(0)\) | \(\mathscr{O}_2(1)\) | \(?\) | \(\mathscr{O}_2(1)\) |
Sean \(\mathscr{O}_i(1)\) y \(\mathscr{O}_i(0)\) los resultados del individuo \(i\) bajo el tratamiento y en ausencia de este, respectivamente. En el modelo tradicional, se asume que no existe interferencia entre individuos, es decir, que la variable de resultado \(\mathscr{O}_i\) depende únicamente de su asignación al grupo de tratamiento o al de control, y es independiente de las asignaciones de las demás unidades experimentales. Así, se define el efecto causal para el individuo \(i\) (también llamado efecto causal a nivel de unidad) como una función únicamente de \(\mathscr{O}_i(1)\) y \(\mathscr{O}_i(0)\)
\[ \tau_{i}=\left[\mathscr{O}_i(1) - \mathscr{O}_i(0)\right] \]
Y a nivel de la población, como \[ \tau_{ATE}=\frac{1}{N}\sum_{i=1}^{N} \left[\mathscr{O}_i(1) - \mathscr{O}_i(0)\right] \]
2.1 Inferencia Causal
La asignación de los individuos al grupo de tratamiento y al de control puede representarse mdiante un vector. Sea \(Z = (Z_1, \ldots, Z_{N})^T \in \{0, 1\}^{N}\) el vector aleatorio de asignación de la intervención, donde \(N\) es el número total de individuos; si la unidad experimental \(i\) recibe el tratamiento, \(Z_i\) es igual a 1, de lo contrario es igual a 0.
Por diseño experimental entenderemos una manera de elegir dicho vector y \(p_z = P(Z = z)\) denotará la probabilidad de que la asignación de tratamiento \(z\) sea generada por el diseño.
Considerando un experimento completamente aleatorizado, en el cual \(T\) individuos son asignados al grupo de tratamiento, y \(C = N - T\) al de control, se tiene que la probabilidad de que la asignación del tratamiento \(z\) sea elegida entre todas las posibles es igual a
\[p_z = \binom{N}{T} ^{-1},\]
Asimismo, se tiene que cada componente \(Z_i\) del vector \(Z\) es una variable aleatoria sigue una distribución Bernoulli con parámetro \(p_{i}=P(Z_i=1)\). Ahora bien, a partir del diseño, podemos determinar la probabilidad de que el individuo \(i\) sea seleccionado en el grupo de tratamiento, \(p_{i}\), de la siguiente manera,
\[p_{i} = \binom{N}{T} ^{-1}\binom{N-1}{T-1}=\frac{T}{N}, \forall i \in V\]
y la probabilidad de que \(i\) y \(j\) lo sean simultáneamente, \(p_{ji}=P(Z_{i}=1, Z_{j}=1)\),
\[p_{ji} = \binom{N}{T} ^{-1}\binom{N-2}{T-2}=\frac{T(T-1)}{N(N-1)}, \forall i \ne j \in V\]
Por tanto, se tiene que
- \(E(Z_i) = p_{i} = \frac{T}{N}\)
- \(\text{Var}(Z_i) = p_{i}(1-p_{i}) = \frac{TC}{N^2}\)
- \(\text{Cov}(Z_i, Z_j) = p_{ij} - p_i p_j = \frac{T(T-1)}{N(N-1)}-\frac{T^2}{N^2}=\frac{-T(N-T)}{N^2(N-1)}, \quad \forall i \ne j \in V\)
Para obtener un estimador insesgado de \(\tau_{ATE}\), ajustamos el total muestral ponderando por el inverso de la probabilidad de inclusión al grupo de tratamiento \(p_{i}\) y
\[ \hat{\tau}_{ATE}(z) = \frac{1}{N} \sum_{i=1}^{N} \left[ \frac{z_i \mathscr{O}_i(1)}{p_{i}} - \frac{(1 - z_i) \mathscr{O}_i(0)}{1-p_{i}} \right] \]
\[ \hat{\tau}_{ATE}(z) = \frac{1}{N} \sum_{i=1}^{N} \left[ \frac{z_i \mathscr{O}_i(1)}{T / N} - \frac{(1 - z_i) \mathscr{O}_i(0)}{C / N} \right] \]
y su varianza por,
\[ V(\hat{\tau}_{ATE}(z)) = V( \frac{1}{N} \sum_{i=1}^{N} \left[ \frac{z_i}{p_{i}}\mathscr{O}_i(1) - \frac{(1 - z_i)}{1-p_{i}} \mathscr{O}_i(0) \right] ) \]
\[ = V(\frac{1}{N} \sum_{i=1}^{N} \left[ (\frac{\mathscr{O}_i(1)}{p_{i}}+ \frac{\mathscr{O}_i(0)}{1-p_{i}})z_i- \frac{\mathscr{O}_i(0)}{1-p_{i}} \right]) \]
\[ =\frac{1}{N^2}\sum_i \left[ \left( \frac{\mathscr{O}_i(1)}{p_{i}} + \frac{\mathscr{O}_i(0)}{1-p_{i}} \right)^2 V(Z_i) \right] + \sum_{i \ne j} \sum \left[ \left( \frac{\mathscr{O}_i(1)}{p_{i}} + \frac{\mathscr{O}_i(0)}{1-p_{i}} \right) \left( \frac{\mathscr{O}_j(1)}{p_{j}} + \frac{\mathscr{O}_j(0)}{1-p_{i}} \right) \text{Cov}(Z_i, Z_j) \right] \]
Finalmente, se obtiene que
\[Var(\hat{\tau}_{ATE}(\mathbf{Z})) = \frac{S_c^2}{C} + \frac{S_t^2}{T} - \frac{S_{tc}^2}{N}\] donde \[S_c^2 = \frac{1}{N - 1} \sum_{i=1}^{N} \left[ \mathscr{O}_i(1) - \bar{\mathscr{O}}(0) \right]^2 \quad \text{,}\] \[S_t^2 = \frac{1}{N - 1} \sum_{i=1}^{N} \left[ \mathscr{O}_i(1) - \bar{\mathscr{O}}(1) \right]^2 \quad \text{y}\] \[S_{tc}^2 = \frac{1}{N - 1} \sum_{i=1}^{N} \left[ \mathscr{O}_i(1) - \mathscr{O}_i(0) - \tau_{ATE}\right]^2 \]
En la práctica no es posible estimar \(S_{tc}^2\), puesto que no observamos simultáneamente \(\mathscr{O}_i(1)\) y \(\mathscr{O}_i(0)\). Bajo el supuesto de efectos unitarios constantes (\(\mathscr{O}_i(1) - \mathscr{O}_i(0) = \tau_{ATE}\)), se puede obtener un estimador insesgado a partir de las varianzas muestrales.
3 Interferencia de Red
En el contexto de las redes, se puede esperar que haya interferencia, es decir, \[\mathscr{O}_i(\mathbf{Z}) \neq \mathscr{O}_i(Z_i)\]
Bajo este nuevo escenario, la respuesta del individuo \(i\) depende de la exposición completa de dicho individuo al vector de asignación de tratamientos \(\mathbf{z}\), en lugar de sólo el producto del tratamiento \(z_i\) al cual fue asignado ese individuo (como sucede en SUTVA).
3.1 Ejemplo: Red de Huelga
La nueva dirección en una planta forestal propuso cambios en el paquete de compensación para los trabajadores. Inicialmente, dos negociadores sindicales explicaron los cambios a los demás, pero no fueron aceptados, resultando en una huelga. Un consultor externo analizó la red de comunicación entre 24 trabajadores relevantes y se explicaron los cambios a dos de ellos, quienes discutieron con sus colegas. En dos días, los trabajadores pidieron reabrir las negociaciones y la huelga se resolvió.
La red de trabajadores consiste en un grafo simple no dirigido y consta de una única componente conectada. La red final tenía 24 vértices y 38 aristas. Una arista entre dos vértices indica que los empleados se comunicaron con una frecuencia mínima suficiente sobre la huelga.
# Librerías
suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(sand)))
# Cargar datos
data(strike)
upgrade_graph(strike)
# Resumen
summary(strike)## IGRAPH 3ec2cc3 U--- 24 38 --
## + attr: names (v/c), race (v/c)
## + edges from 3ec2cc3:
## [1] 1-- 2 1-- 3 1-- 4 2-- 3 2-- 4 3-- 4 3-- 5 5-- 6 5-- 7 5-- 8
## [11] 5--11 5--12 5--15 6-- 7 7-- 9 8-- 9 8--11 9--10 9--11 11--12
## [21] 11--13 12--13 13--14 14--15 15--16 15--20 15--21 15--24 16--17 17--18
## [31] 17--24 18--19 18--24 19--20 21--22 21--23 21--24 22--23
## IGRAPH 3ec2cc3 U--- 24 38 --
## + attr: names (v/c), race (v/c)Hay tres subgrupos presentes en la red, codificados por el atributo nodal \(\verb|race|\):
| \(\verb|race|\) | Decodificación |
|---|---|
| OE | Empleados mayores de habla inglesa |
| YE | Empleados jóvenes de habla inglesa |
| YS | Empleados jóvenes de habla hispana |
# O -> Older
# Y -> Younger
# E -> English-speaking
# S -> Spanish-speaking
table(V(strike)$race)##
## OE YE YS
## 11 9 4El experimento consiste en explicar los cambios en las condiciones de pago para influir en el comportamiento organizacional, es decir, aprovechando la estructura de comunicación entre un grupo de empleados en huelga se busca conciliar la disputa en una instalación de fabricación de productos forestales.
En la siguiente figura se presentan las relaciones al interior del sindicato, según características nodales:
La pertenencia étnica y el ciclo de vida se indican según la forma de los vértices. El círculo refiere a jóvenes hispanos, cuadrado a jóvenes anglohablantes y triángulo a adultos anglohablantes.
Los individuos seleccionados para la intervención son coloreados en rojo.
Los dos representantes sindicales, negociadores en el intento de comunicación inicial, están marcados con un asterisco (es decir, Sam y Wendle).
# Darle forma de triángulo a los vértices
mytriangle <- function(coords, v=NULL, params) {
vertex.color <- params("vertex", "color")
if (length(vertex.color) != 1 && !is.null(v)) {
vertex.color <- vertex.color[v]
}
vertex.size <- 1/200 * params("vertex", "size")
if (length(vertex.size) != 1 && !is.null(v)) {
vertex.size <- vertex.size[v]
}
symbols(x=coords[,1], y=coords[,2], bg=vertex.color,
stars=cbind(vertex.size, vertex.size, vertex.size),
add=TRUE, inches=FALSE)
}
add_shape("triangle", clip=shapes("circle")$clip,
plot=mytriangle)
# Indicar "race" a través de la forma del vértice
V(strike)[V(strike)$race=="YS"]$shape <- "circle"
V(strike)[V(strike)$race=="YE"]$shape <- "square"
V(strike)[V(strike)$race=="OE"]$shape <- "triangle"
# Distinguir los 4 representantes que difundieron la estrategia mediante colores
nv <- vcount(strike)
z <- numeric(nv)
z[c(5,15,21,22)] <- 1
V(strike)$color <- rep("white",nv)
V(strike)[V(strike)$race=="YS"]$color <- "#19D339"
V(strike)[V(strike)$race=="YE"]$color <- adjustcolor("royalblue", 0.7)
V(strike)[V(strike)$race=="OE"]$color <- "navy"
V(strike)[z==1]$color <- "red3"#Visualización
par(mar = c(0.5, 1, 1, 0.5), mgp = c(1.2, 0.3, 0), mfrow = c(1, 1))
set.seed(42)
my.dist <- c(rep(1.8,4),rep(2.2,9),rep(2,11))
l <- layout_with_kk(strike) # Kamada-Kawai
plot(strike,layout=l,vertex.label=V(strike)$names,
vertex.label.degree=-pi/3,
vertex.label.dist=my.dist,
vertex.label.color="black")
legend("topleft",
legend = c(paste("Orden = ", vcount(strike), sep=""),
paste("Tamaño = ", ecount(strike), sep=""),
paste("Densidad = ", round(edge_density(strike), 2), sep=""),
paste("Grado promedio = ", round(mean(degree(strike)), 2), sep=""),
paste("Diámetro = ", round(diameter(strike, directed = FALSE, unconnected = TRUE), 2), sep=""),
paste("Distancia promedio = ", round(mean_distance(strike, directed = FALSE, unconnected = TRUE), 2), sep=""),
paste("Triángulos = ", sum(count_triangles(strike))/2, sep=""),
paste("Transitividad = ", round(transitivity(strike, type="global"), 2), sep=""),
paste("Num. clique = ", clique_num(graph = strike), sep="")),
bty = "n", cex = 0.7, text.width = 2)
legend("bottomright",
legend = c("adulto, anglohablante",
"joven, anglohablante",
"joven, hispanohablante",
"tratados"),
col = c("navy", adjustcolor("royalblue", 0.7), "#19D339", "red3"),
pch = 16,
pt.cex = 1.5,
bty = "n",
title = "",
cex = 0.7)
En primer lugar, se pueden observar las siguientes características sobre la red:
Los 38 enlaces observados equivalen al 14% de las aristas posibles (densidad = 0.14).
En promedio, cada sindicalista se relaciona con 3.17 trabajadores adicionales.
El 88% de los nodos tienen 4 o menos enlaces.
En promedio, se necesita recorrer un camino de casi 3 pasos para ir de un nodo a otro en la red (distancia media = 2.99).
El par de nodos más alejados entre sí se encuentran a 6 aristas de distancia (diámetro = 6).
# k-conectividad
print("k-conectividad")
vertex_connectivity(strike)
# Puntos de articulación
print("Puntos de articulación")
V(strike)[articulation_points(strike)]$names## [1] "k-conectividad"
## [1] 1
## [1] "Puntos de articulación"
## [1] "Sam*" "Norm" "Gill" "Bob" "Alejandro"Al analizar la conectividad del grafo, se puede constatar que basta con remover un solo vértice bien elegido para dividir la componente gigante en subgrafos más pequeños. Por lo tanto, decimos que la red del sindicato es 1-conectada. Esto implica que el flujo de información entre los trabajadores puede ser interrumpido al eliminar al menos uno de los nodos de articulación de la red.
centrality_table <- data.frame(
`x1` = V(strike)[order(degree(strike, normalized = TRUE), decreasing = TRUE)[1:5]]$names,
`x2`= V(strike)[order(closeness(strike, normalized = TRUE), decreasing = TRUE)[1:5]]$names,
`x3` = V(strike)[order(betweenness(strike, normalized = TRUE), decreasing = TRUE)[1:5]]$names,
`x4` = V(strike)[order(eigen_centrality(strike, directed = FALSE, scale = TRUE)$vector, decreasing = TRUE)[1:5]]$names
)
colnames(centrality_table) <- c( "Grado", "Cercanía", "Intermediación", "Vector propio")
# Mostrar la tabla
kable(centrality_table, caption = "Nodos destacados, grafo Strike") %>%
kable_styling(full_width = FALSE, position = "center") %>%
add_header_above(c("Medidas de Centralidad, Top 5" = 4))| Grado | Cercanía | Intermediación | Vector propio |
|---|---|---|---|
| Bob | Bob | Bob | Bob |
| Norm | Norm | Norm | John |
| John | John | Alejandro | Hal |
| Alejandro | Ultrecht | Sam* | Lanny |
| Gill | Alejandro | Ultrecht | Norm |
print("betweenness")
V(strike)[z==1]$names
rank(-betweenness(strike))[z==1]## [1] "betweenness"
## [1] "Bob" "Norm" "Sam*" "Wendle*"
## [1] 1 2 4 21print("closeness")
V(strike)[z==1]$names
rank(-closeness(strike))[z==1]## [1] "closeness"
## [1] "Bob" "Norm" "Sam*" "Wendle*"
## [1] 1.0 2.0 6.5 22.0Bajo las métricas estándar, los individuos más centrales de la red son Norm, Bob y John, a la vez que Frank, Xavier y Wendle* son los individuos más periféricos. La centralidad de Bob y Norm justifica que hayan sido elegidos luego de analizar la estructura de la red, para propagar la información hacia más sectores. Al contrario, Sam y Wendle, los representantes sindicales, parecen alejados del resto de trabajadores.
# Comunidades- fast_greedy
par(mar = c(0.5, 1, 1, 0.5), mgp = c(1.2, 0.3, 0), mfrow = c(1, 1))
set.seed(42)
kc <- cluster_fast_greedy(strike)
my.dist <- c(rep(1.8,4),rep(2.2,9),rep(2,11))
l <- layout_with_kk(strike) # Kamada-Kawai
plot(kc,strike,layout=l,vertex.label=V(strike)$names,
vertex.label.degree=-pi/3,
vertex.label.dist=my.dist,
vertex.label.color="black")
legend("bottomright", legend = unique(V(strike)$race), pt.bg = "white",
pch = c(21, 22, 24), bty = "n", col = "#777777", pt.cex = 1.5, pt.lwd = 2)
print("Nodos en el clique máximo")
V(strike)[maximal.cliques(strike)[[which.max(sapply(maximal.cliques(strike), length))]]]$names## Warning: `maximal.cliques()` was deprecated in igraph 2.0.0.
## ℹ Please use `max_cliques()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## [1] "Nodos en el clique máximo"
## [1] "Alejandro" "Domingo" "Eduardo" "Carlos"En relación con la cohesión del grafo, se tiene que:
Aunque no existen regiones aisladas en la red, las regiones de máxima conectividad son pequeñas.
El subgrupo más altamente conectado es relativamente pequeño, siendo el clique máximo de 4 vértices.
El clique máximo corresponde a la comunidad de habla hispana.
Al utilizar tres algortimos para la segmentación del grafo: agrupación jerárquica (fast-greedy), particionamiento espectral y el algoritmo de Girvan-Newman (edge-betweenness), se evidenció que las cuatro comunidades detectadas están claramente delimitadas, siendo así que:
Los dos subgrupos más grandes reúnen al 70% de los nodos, con 10 y 7 vértices respectivamente.
Los dos subgrupos más pequeños están formados por 4 o menos individuos cada uno.
4 Inferencia Causal bajo Interferencia de Red
Al buscar extender el marco clásico de resultados potenciales a experimentos en red, surge una tensión entre la complejidad de la interferencia y su impacto en:
Definir parámetros.
Producir estimadores y cuantificar la incertidumbre asociada.
4.1 Modelos de Exposición en Redes
Ahora la respuesta del individuo \(i\) se entiende como el resultado de la exposición completa de dicho individuo al vector de asignación de tratamientos \(\mathbf{z}\), en lugar de sólo el producto del tratamiento \(z_i\) al cual fue asignado ese individuo (como sucede en SUTVA).
El problema es que habrá \(2^{N}\) posibles exposiciones por cada uno de los \(N\) individuos, haciendo imposible el inferencia causal.
4.1.1 Mapeos de Exposición
Es una restricción de modelado, propuesta por Aronow y Samii, sobre el grado en que la interferencia de otros individuos en la red afecta la exposición de un individuo \(i\).
Se asume que efectivamente hay un número finito \(K\) de condiciones \(\{c_1, \ldots, c_K\}\) al que el individuo \(i\) está expuesto.
Decimos que \(i\) está expuesto a la condición \(k\) si \[f(\mathbf{z}, \mathbf{x}_i) = c_k\]
donde
\(\mathbf{z}\) es el vector de asignación de tratamientos.
\(\mathbf{x}_i\) es el vector de la i-esíma columna de la matriz de adyacencia.
\(f\) es la función del mapeo de exposición.
Denotamos \(A\) como la matriz de adyacencia de la red \(G\), y el vector \(\mathbf{x}_i\) como la i-ésima columna de esta matriz. Entonces, definimos la categorización de exposición simple de cuatro niveles: \[f(\mathbf{z}, \mathbf{x}_i) = \begin{cases} c_{11}(\text{Exposición Directa + Indirecta}), & z_i\mathit{I}_{ \{\mathbf{z}^T\mathbf{x}_i>0\} } = 1\\ c_{10}(\text{Exposición Directa Aislada}), & z_i\mathit{I}_{ \{\mathbf{z}^T\mathbf{x}_i=0\} } = 1\\ c_{01}(\text{Exposición Indirecta}), & (1 - z_i)\mathit{I}_{ \{\mathbf{z}^T\mathbf{x}_i>0\} } = 1\\ c_{00}(\text{No Exposición}), & (1 - z_i)\mathit{I}_{ \{\mathbf{z}^T\mathbf{x}_i=0\} } = 1 \end{cases}\] donde \(\mathbf{z}^T\mathbf{x}_i\) es igual al número de vecinos directamente expuestos del individuo \(i\).
4.1.2 Efectos del Tratamiento bajo Interferencia de Red
La adición de la interferencia al marco de resultados potenciales requiere refinar la noción de efectos causales. El efecto promedio del tratamiento (ATE) bajo interferencia de la red se generaliza como: \[ \tau_{ATE} = \frac{1}{N} \sum_{i=1}^{N} \left[ \mathscr{O}_i(\mathbf{1}) - \mathscr{O}_i(\mathbf{0}) \right] = \bar{\mathscr{O}}(\mathbf{1}) - \bar{\mathscr{O}}(\mathbf{0}) \] donde \(\mathbf{1}\) y \(\mathbf{0}\) son ahora vectores de longitud \(N\) de unos y ceros, respectivamente. Compara resultados potenciales bajo tratamiento completo (\(\mathbf{z=1}\)) versus control total (\(\mathbf{z=0}\)).
El problema es que este parámetro representa una agregación de ambos efectos (directos e indirectos) de la asignación del tratamiento y, en contextos donde se buscan efectos “puros” del tratamiento, la interferencia puede ser una molestia.
Bajo el marco del mapeo de exposición de Aronow y Samii, usamos la diferencia \(\mathscr{O}_i(c_k) - \mathscr{O}_i(c_l)\) para representar el efecto causal de la condición de exposición \(k\) versus la \(l\) para el individuo \(i\). Entonces, el efecto causal promedio de exposición a la condición \(k\) vesus \(l\) se define como: \[ \tau(c_k, c_l) = \frac{1}{N} \sum_{i=1}^{N} \left[ \mathscr{O}_i(c_k) - \mathscr{O}_i(c_l) \right] = \bar{\mathscr{O}}(c_k) - \bar{\mathscr{O}}(c_l) \]
Entre los niveles de exposición posibles: \(c_{11}, c_{10}, c_{01}, c_{00}\), podemos definir los siguientes contrastes:
\(\tau(c_{01}, c_{00})\) que captura el efecto general indirecto del tratamiento.
\(\tau(c_{10}, c_{00})\) que captura el efecto general directo del tratamiento.
\(\tau(c_{11}, c_{00})\) que captura el efecto total del tratamiento.
4.1.3 Ejemplo: Red de Huelga
Para poder identificar la exposición de cada individuo bajo el modelo de exposición propuesto por Aronow y Samii, se debe calcular el indicador de vecinos directamente expuestos:
A <- as_adjacency_matrix(strike)
I.ex.nbrs <- as.numeric(z%*%A > 0)
I.ex.nbrs## [1] 0 0 1 0 1 1 1 1 0 0 1 1 0 1 1 1 0 0 0 1 1 1 1 1En la red de huelga, hay cuatro empleados que recibieron exposición tanto directa como indirecta \(c_{11}\),
V(strike)[z*I.ex.nbrs==1]$exposure <- "c[11]"
V(strike)[z*I.ex.nbrs==1]$names## [1] "Bob" "Norm" "Sam*" "Wendle*"once empleados recibieron sólo exposición indirecta \(c_{01}\)
V(strike)[(1-z)*I.ex.nbrs==1]$exposure <- "c['01']"
V(strike)[(1-z)*I.ex.nbrs==1]$names## [1] "Alejandro" "Mike" "Ike" "Hal" "John" "Lanny"
## [7] "Ozzie" "Paul" "Vern" "Xavier" "Ultrecht"y ninguno recibió exposición directa aislada \(c_{10}\).
V(strike)[z*(1-I.ex.nbrs)==1]$exposure <- "c[10]"
V(strike)[z*(1-I.ex.nbrs)==1]$names## character(0)El resto de los empleados no fueron expuestos \(c_{00}\).
V(strike)[(1-z)*(1-I.ex.nbrs)==1]$exposure <- "c['00']"
V(strike)[(1-z)*(1-I.ex.nbrs)==1]$names## [1] "Domingo" "Carlos" "Eduardo" "Gill" "Frank" "Karl" "Quint"
## [8] "Russ" "Ted"table(V(strike)$exposure)##
## c['00'] c['01'] c[11]
## 9 11 4La categorización de exposición simple de cuatro niveles se puede visualizar como
suppressMessages(suppressWarnings(library(RColorBrewer)))
# Definir colores para los valores de exposición
set.seed(42)
palette_colors <- brewer.pal(length(unique(V(strike)$exposure)), "Set3")
exposure_colors <- setNames(palette_colors, unique(V(strike)$exposure))
vertex_colors <- exposure_colors[V(strike)$exposure]
# Visualización
plot(strike,layout = l, vertex.label = V(strike)$names,
vertex.label.degree = -pi/3, vertex.label.dist = my.dist,
vertex.color = vertex_colors)
legend("topright", legend = parse(text = unique(V(strike)$exposure)),
pt.bg = exposure_colors, bty = "n", pch = 21, col = "#777777")
Supongamos que la reacción de los empleados se califica en una escala de 1 (no receptivo) a 10 (completamente receptivo), registrando que la receptividad incrementa según la exposición. Para cada individuo \(i \in V\), se tiene que:
\(\mathscr{O}_i(c_{00}) = 1\): mínima receptividad sin exposición.
\(\mathscr{O}_i(c_{01}) = 5\): receptividad aumenta con solo exposición indirecta.
\(\mathscr{O}_i(c_{10}) = 7\): mayor receptividad con exposición directa aislada.
\(\mathscr{O}_i(c_{11}) = 10\): máxima receptividad con exposición completa.
Por tanto, los correspondientes efectos causales son \(\tau(c_{01}, c_{00}), \tau(c_{10}, c_{00}), \tau(c_{11}, c_{00})\), que toman respectivamente los siguientes valores.
O.c00 <- 1.0
O.c01 <- 5
O.c10 <- 7
O.c11 <- 10.0
c(O.c01, O.c10, O.c11) - O.c00## [1] 4 6 94.2 Efectos Causales bajo Interferencia de Red
La elección de la asignación de tratamientos es un elemento clave del diseño experimental de redes. En un modelo clásico de inferencia causal, la asignación se captura a través de la distribución \(p_{\mathbf{z}}\), es decir, la probabilidad de asignación del tratamiento \(\mathbf{z}\).
Bajo interferencia, es más adecuado considerar la exposición general de los individuos al tratamiento. Por tanto, para caracterizar un diseño experimental de red, es útil cuantificar cómo la asignación de tratamiento induce exposición bajo un modelo específico de exposición en la red.
En el contexto de los mapeos de exposición en redes, la probabilidad de un individuo \(i\) sea sujeto a la condición de exposición \(k\) se define como \[p_{i}^{e}(c_k) = \sum_{\mathbf{z}} p_{\mathbf{z}} \mathit{I}_{ \{ f(\mathbf{z}, \mathbf{x}_i) = c_k \} }\] Para todos los individuos \(i\) y pares de individuos \(i\) y \(j\), los valores de las distintas probabilidades de exposición pueden recuperarse a partir de dos clases de matrices.
Supongamos que hay \(\mathit{M}\) posibles asignaciones de tratamiento \(\mathbf{z}\), siendo \[\begin{equation} \mathit{I}_{k} = \begin{bmatrix} \mathit{I}_{\{f(\mathbf{z}_1, \mathbf{x}_1) = c_k\}} & \mathit{I}_{\{f(\mathbf{z}_2, \mathbf{x}_1) = c_k\}} & \cdots & \mathit{I}_{\{f(\mathbf{z}_{\mathit{M}}, \mathbf{x}_1) = c_k\}} \\ \mathit{I}_{\{f(\mathbf{z}_1, \mathbf{x}_2) = c_k\}} & \mathit{I}_{\{f(\mathbf{z}_2, \mathbf{x}_2) = c_k\}} & \cdots & \mathit{I}_{\{f(\mathbf{z}_{\mathit{M}}, \mathbf{x}_2) = c_k\}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathit{I}_{\{f(\mathbf{z}_1, \mathbf{x}_{N}) = c_k\}} & \mathit{I}_{\{f(\mathbf{z}_2, \mathbf{x}_{N}) = c_k\}} & \cdots & \mathit{I}_{\{f(\mathbf{z}_{\mathit{M}}, \mathbf{x}_{N}) = c_k\}} \end{bmatrix} \end{equation}\]
Además, sea \(\mathbf{P} = diag(p_{z_1},\ldots, p_{z_\mathit{M}})\). Entonces, en la primera matriz simétrica \(N \times N\), podemos recuperar para una condición de exposición fija \(c_k\), tanto la probabilidad de exposición individual \(p_{i}^{e}(c_k)\) como las probabilidades de exposición conjunta \(p_{ij}^{e}(c_k)\), para todos los individuos \(i, j = 1,\ldots,N\): \[\begin{equation} \mathit{I}_k \mathbf{P} \mathit{I}_{k}^{T} = \begin{bmatrix} p_{1}^{e}(c_k) & p_{12}^{e}(c_k) & \cdots & p_{1N}^{e}(c_k) \\ p_{21}^{e}(c_k) & p_{2}^{e}(c_k) & \cdots & p_{2N}^{e}(c_k) \\ \vdots & \vdots & \ddots & \vdots \\ p_{N1}^{e}(c_k) & p_{N2}^{e}(c_k) & \cdots & p_{N}^{e}(c_k) \end{bmatrix} \end{equation}\] y la segunda matriz no simétrica \(N \times N\) nos permite recuperar las probabilidades de exposición conjunta \(p_{ij}^{e}(c_k, c_l)\) para todos los pares de individuos \(i\) y \(j\), para condiciones de exposición fijas \(c_k\) y \(c_l\): \[\begin{equation} \mathit{I}_k \mathbf{P} \mathit{I}_{l}^{T} = \begin{bmatrix} 0 & p_{12}^{e}(c_k, c_l) & \cdots & p_{1N}^{e}(c_k, c_l) \\ p_{21}^{e}(c_k, c_l) & 0 & \cdots & p_{2N}^{e}(c_k, c_l) \\ \vdots & \vdots & \ddots & \vdots \\ p_{N1}^{e}(c_k, c_l) & p_{N2}^{e}(c_k, c_l) & \cdots & 0 \end{bmatrix} \end{equation}\] Note que la diagonal está compuesta por ceros debido al supuesto de que los individuos solo pueden caer en una categoría de exposición.
El problema es que en la práctica, es probable que las dos matrices sean computacionalmente inviables para calcular directamente \(p_{i}^{e}(c_k)\), \(p_{ij}^{e}(c_k)\) y \(p_{ij}^{e}(c_k, c_l)\).
En el marco del mapeo de exposición, las estimaciones de interés son:
\(\bar{\mathscr{O}}(c_k)\) siendo la resultado promedio potencial bajo la condición de exposición \(c_k\) para \(k = 1,\ldots,K\).
\(\tau(c_k, c_l) = \bar{\mathscr{O}}(c_k) - \bar{\mathscr{O}}(c_l)\) siendo el efecto causal promedio de la condición de exposición \(c_k\) versus \(c_l\).
Note que los \(\bar{\mathscr{O}}(c_k)\) corresponderán a una muestra de probabilidades desiguales sin reemplazo de resultados potenciales bajo una condición dada, y las estimaciones del 2 son todas funciones lineales simples de las del 1.
Suponiendo que conocemos o podemos aproximarnos con precisión a las probabilidades de exposición relevantes, la inferencia se realiza siguiendo el método de Horvitz y Thompson.
Este estimador insesgado y bien definido para \(\bar{\mathscr{O}}(c_k)\), considera el muestro de probabilidades desiguales mediante el uso de ponderación de probabilidad inversa.
\[\hat{\bar{\mathscr{O}}}(c_k) = \frac{1}{N} \sum_{i=1}^{N} \mathit{I}_{ \{ f(\mathbf{z}, \mathbf{x}_i) = c_k \} } \frac{\mathscr{O}_i(c_k)}{p_i^e(c_k)}\]
A su vez, \(\hat{\tau}(c_k, c_l) = \hat{\bar{\mathscr{O}}}(c_k) - \hat{\bar{\mathscr{O}}}(c_l)\) es un estimador insesgado para \(\tau(c_k, c_l)\).
Las varianzas de estos estimadores son, respectivamente,
\[Var[\hat{\bar{\mathscr{O}}}(c_k)] = \frac{1}{N^2} \Biggl\{ \sum_{i=1}^{N} p_i^e(c_k) [1-p_i^e(c_k)] \biggl[ \frac{\mathscr{O}_i(c_k)}{p_i^e(c_k)} \biggl]^2 + \sum_{i=1}^{N}\sum_{i\neq1} [p_{ij}^e(c_k) - p_i^e(c_k)p_j^e(c_k)] \frac{\mathscr{O}_i(c_k)}{p_i^e(c_k)} \frac{\mathscr{O}_j(c_k)}{p_j^e(c_k)} \Biggl\}\] \[Var(\hat{\tau}(c_k, c_l)) = Var[\hat{\bar{\mathscr{O}}}(c_k)] + Var[\hat{\bar{\mathscr{O}}}(c_l)] - 2Cov[\hat{\bar{\mathscr{O}}}(c_k), \hat{\bar{\mathscr{O}}}(c_l)]\]
4.2.1 Ejemplo: Red de Huelga
Se pueden usar simulaciones de Monte Carlo para aproximar el valor de las probabilidades con precisión arbitraria.
Simulando \(n\) valores de \(p_{\mathbf{z}}\), para cada una de las \(K\) condiciones de exposición podemos entonces formar \(K\) matrices \(\hat{\mathit{I}}_k\) de dimensión \(N \times n\).
Supongamos que en la red de huelga, se aplicó “tratamiento” a cuatro empleados elegidos uniformemente al azar para que desempeñen el rol de negociadores.
# Inicializar
set.seed(41)
m <- 4 # Número de representantes
n <- 10000 # Número de ensayos de Montecarlo
# 4 condiciones de exposición
I11 <- matrix(,nrow=nv,ncol=n)
I10 <- matrix(,nrow=nv,ncol=n)
I01 <- matrix(,nrow=nv,ncol=n)
I00 <- matrix(,nrow=nv,ncol=n)
# Muestreo Monte Carlo
for(i in 1:n){
z <- rep(0,nv) # Vector z
reps.ind <- sample((1:nv), m, replace = FALSE) # Seleccionar representantes
z[reps.ind] <- 1 # Asignar 1 al vector z
reps.nbrs <- as.numeric(z%*%A > 0) # Indicador de vecinos expuestos
# Una matriz por cada condición de exposición
I11[,i] <- z*reps.nbrs
I10[,i] <- z*(1-reps.nbrs)
I01[,i] <- (1-z)*reps.nbrs
I00[,i] <- (1-z)*(1-reps.nbrs)
}Los estimadores insesgados para las dos matrices con las probabilidades de exposición son \(\hat{\mathit{I}}_k\hat{\mathit{I}}_k^T/n\) y \(\hat{\mathit{I}}_k\hat{\mathit{I}}_l^T/n\), y convergen con seguridad por la ley de los grandes números.
I11.11 <- I11%*%t(I11)/n
I10.10 <- I10%*%t(I10)/n
I01.01 <- I01%*%t(I01)/n
I00.00 <- I00%*%t(I00)/nsuppressMessages(suppressWarnings(library(fields)))
# Visualización
names.w.space <- paste(V(strike)$names," ",sep = "")
my.cex.x <- 0.75
my.cex.y <- 0.75
plot_Ik <- function(M, title) {
# Márgenes
par(mar = c(4,4,4,4))
# Matriz
image(M, zlim = c(0, 0.7), xaxt = "n", yaxt = "n", col = cm.colors(16))
# Nombres ejes
mtext(side = 1, text = names.w.space, at = seq(0.0, 1.0, (1/23)),
las = 3, cex = my.cex.x)
mtext(side = 2, text = names.w.space, at = seq(0.0, 1.0, 1/23),
las = 1, cex = my.cex.y)
# Título
mtext(side = 3, text = title, at = 0.5, las = 1)
# Agregar líneas para diferenciar grupos
u <- 1/23
uo2 <- 1/46
xmat <- cbind(rep(3*u+uo2,2), rep(12*u+uo2,2))
ymat <- cbind(c(0-uo2,1+uo2), c(0-uo2,1+uo2))
matlines(xmat, ymat, lty = 1, lw = 1, col = "black")
matlines(ymat, xmat, lty = 1, lw = 1, col = "black")
# Añadir barra de colores
image.plot(legend.only = TRUE, zlim = range(seq(0.0,0.7,0.1)), col = cm.colors(16))
}
# Crear los cuatro gráficos
par(mfrow = c(2, 2))
plot_Ik(I11.11, expression("Exposición Directa + Indirecta"~(c["11"])))
plot_Ik(I10.10, expression("Exposición Directa Aislada"~(c["10"])))
plot_Ik(I01.01, expression("Exposición Indirecta"~(c["01"])))
plot_Ik(I00.00, expression("No Exposición"~(c["00"])))
Con solo cuatro de 24 actores tratados como negociadores, solo las
probabilidades de exposición indirecta y sin exposición tienen valores
altos. Se observa que en el subgrupo de empleados jóvenes de habla
hispana hay alta probabilidad de no exposición.
Al suponer que los cuatro representantes de la red de huelga fueron elegidos al azar, las estimaciones de los efectos causales promedio \(\tau(c_{11}, c_{00})\), \(\tau(c_{10}, c_{00})\) y \(\tau(c_{01}, c_{00})\) se obtienen a continuación.
# Vector de asignación del tratamiento
z <- rep(0,nv)
z[c(5,15,21,22)] <- 1
# Indicador de vecinos expuestos
reps.nbrs <- as.numeric(z%*%A > 0)
# Condiciones de exposición
c11 <- z*reps.nbrs
c10 <- z*(1-reps.nbrs)
c01 <- (1-z)*reps.nbrs
c00 <- (1-z)*(1-reps.nbrs)
# Estimadores Horvitz y Thompson
Obar.c11 <- O.c11*mean(c11/diag(I11.11))
Obar.c10 <- O.c10*mean(c10/diag(I10.10))
Obar.c01 <- O.c01*mean(c01/diag(I01.01))
Obar.c00 <- O.c00*mean(c00/diag(I00.00))
print(c(Obar.c11, Obar.c10, Obar.c01) - Obar.c00)## [1] 21.4173702 -0.8001775 5.9828283Estos valores se comparan bastante mal con sus respectivos valores objetivo 9, 6 y 4, sugiriendo que pese a ser insesgado, la varianza de estos estimadores sea motivo de preocupación.
El problema es que \(Var(\hat{\tau}(c_k, c_l))\) no se puede estimar de manera insesgada o consistente, debido a su último término que depende de resultados potenciales (no observados).
Es posible usar la simulación de Monte Carlo para producir un estimador de esta varianza con sesgo conservador y obtener una idea del rendimiento de estos estimadores.
set.seed(41)
n <- 10000
Obar.c11 <- numeric()
Obar.c10 <- numeric()
Obar.c01 <- numeric()
Obar.c00 <- numeric()
for(i in 1:n){
z <- rep(0,nv)
reps.ind <- sample((1:nv),m,replace=FALSE)
z[reps.ind] <- 1
reps.nbrs <- as.numeric(z%*%A > 0)
c11 <- z*reps.nbrs
c10 <- z*(1-reps.nbrs)
c01 <- (1-z)*reps.nbrs
c00 <- (1-z)*(1-reps.nbrs)
Obar.c11 <- c(Obar.c11, O.c11*mean(c11/diag(I11.11)))
Obar.c10 <- c(Obar.c10, O.c10*mean(c10/diag(I10.10)))
Obar.c01 <- c(Obar.c01, O.c01*mean(c01/diag(I01.01)))
Obar.c00 <- c(Obar.c00, O.c00*mean(c00/diag(I00.00)))
}
# Average causal effects
ACE <- list(Obar.c11-Obar.c00,
Obar.c10-Obar.c00,
Obar.c01-Obar.c00)
print(sapply(ACE, mean))## [1] 9 6 4# Diferencias con el valor objetivo
print(sapply(ACE, mean) - c(O.c11-O.c00, O.c10-O.c00, O.c01-O.c00))## [1] 0 0 0# Errores estándar
print(sapply(ACE, sd))## [1] 8.928154 3.701966 1.590345# Coeficientes de variación
sapply(ACE, sd) / c(9, 6, 4)## [1] 0.9920171 0.6169943 0.3975862Vemos que la diferencia entre el valor esperado y los estimadores de los tres efectos causales promedio están cercanos a 0, sin embargo, los errores estándar correspondientes son bastante dispares. En este orden, los estimadores son mejores para capturar \(\tau(c_{01}, c_{00})\), luego \(\tau(c_{10}, c_{00})\) y, por último, \(\tau(c_{11}, c_{00})\).
En conclusión, el diseño experimental permite observar no solo el efecto directo de la intervención sobre los individuos tratados, sino también los efectos indirectos en sus colegas con los que están conectados. Este enfoque es particularmente útil en situaciones donde las decisiones y comportamientos de un individuo pueden impactar significativamente en el comportamiento de otros a través de su red de comunicación. Cambiar de un diseño experimental a un diseño no experimental puede afectar significativamente la capacidad para establecer relaciones causales y generalizar los resultados.
5 Ejemplo 2: Climas Cambiantes de Conflicto
La colección de datos proviene de un experimento de campo a gran escala en 56 escuelas secundarias públicas de Nueva Jersey con 24.191 estudiantes. Se investigó cómo el comportamiento de los individuos influye en sus acciones al observar a personas de su comunidad. La intervención anticconflicto se aplicó aleatoriamente a grupos de 20 a 32 estudiantes en las escuelas seleccionadas. Los datos incluyen redes sociales de los estudiantes, encuestas sobre normas de conflicto y variables demográficas como género, edad, raza, y características del hogar, con un total de 482 variables para 24.471 casos. Las encuestas se realizaron al inicio y al final del año escolar 2012-2013.
Nota sobre acceso a datos: La recopilación de datos no pertenece a las autoras del presente cuaderno, sino que está asociada con la publicación “Changing climates of conflict: A social network experiment in 56 schools” por Elizabeth Levy Paluck, Hana Shepherd, y Peter M. Aronow. Los usuarios pueden visitar el sitio web de la publicación y consultar los paquetes (.zip) incluidos en la documentación de este estudio. Así mismo, la versión de los datos utilizados en este ejercicio se puede obtener por medio del siguiente enlace: https://www.icpsr.umich.edu/web/civicleads/studies/37070.
# Librerías
suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(RColorBrewer)))
suppressMessages(suppressWarnings(library(fields)))
# Cargar datos
load("./Data/DS0001/37070-0001-Data.rda")
# Renombrar datos
assign("conflict", get("da37070.0001"))
rm(da37070.0001)Para el presente ejercicio, con el objetivo de delimitar el análisis, se selecciona al azar una escuela que haya sido tratada. Dentro de esta escuela, se encuentran 3 subgrupos presentes en la red:
- (0) Sin tratamiento ni control: Estudiantes que no recibieron tratamiento.
- (1) Tratamiento: Estudiantes que participaron en la intervención anticonflicto.
- (2) Control: Estudiantes sin tratamiento, ni compañeros tratados que sean referentes sociales.
Para analizar las relaciones entre los estudiantes, se les planteó la siguiente pregunta: “En las últimas semanas, decidí pasar tiempo con estos estudiantes de mi escuela (ya sea en la escuela, fuera de ella o en línea)”. Cada estudiante debía seleccionar un máximo de 5 compañeros.
Asimismo, se utilizaron los eventos disciplinarios reportados por los docentes para cada estudiante como medida de la variable de conflictos; dicha variable fue extraída de los registros administrativos correspondientes a todo el año.
# Escoger al azar una escuela que haya recibido tratamiento
set.seed(123)
num <- sample(unique(conflict$SCHID[conflict$SCHTREAT == "(1) Treatment school (Roots meetings 2012-2013)"]), 1)
# Filtrar por escuela seleccionada
conflict_29 <- conflict[conflict$SCHID == num,]
# Modificar df para crear al grafo
conflict_29 <- conflict_29 %>%
pivot_longer(cols = ST1:ST5, names_to = "ST_names", values_to = "to") %>%
dplyr::mutate(from = as.numeric(ID),
event_disc = rowSums(select(., matches("^DE[0-9]")))) %>%
dplyr::select(c(from, to, TREAT, event_disc)) %>% na.omit()
# Extraer vertices y aristas
vertices <- conflict_29[c(1,3:4)] %>% distinct()
conflict_29 <- conflict_29 %>% filter(to %in% vertices$from)
edges <- conflict_29[1:2]
vertices %>%
group_by(TREAT) %>%
summarise(n = n_distinct(from))## # A tibble: 3 × 2
## TREAT n
## <fct> <int>
## 1 (0) Not treatment or control 217
## 2 (1) Treatment 24
## 3 (2) Control 21En el caso de la escuela 29, la muestra está compuesta por 262 estudiantes, distribuidos de la siguiente forma: 217 estudiantes pertenecen al grupo de sin tratamiento ni control, 24 al grupo de tratamiento, y 21 al grupo de control.
En la siguiente figura se visualiza la estructura de la red, que cuenta con un total de 863 enlaces entre los nodos. La red presenta una densidad de 0.03, lo que indica que solo el 3% de todas las posibles conexiones están presentes. Así mismo, cada nodo está conectado a aproximadamente 7 otros nodos. Por último, en promedio, cualquier par de nodos está separado por poco más de 3 enlaces.
# Grafo
g <- graph_from_data_frame(edges, directed = FALSE, vertices = vertices)
# Layout
set.seed(123)
l <- layout.fruchterman.reingold(g)
# Visualización
palette_colors <- brewer.pal(length(unique(V(g)$TREAT)), "Set2")
exposure_colors <- setNames(palette_colors, unique(V(g)$TREAT))
vertex_colors <- exposure_colors[V(g)$TREAT]
plot(g, layout = l, vertex.label = NA, vertex.color = vertex_colors, vertex.size = 4)
legend("topleft",
legend = names(exposure_colors),
fill = exposure_colors,
title = "TREAT")
legend("bottomleft",
legend = c(paste0("Orden = ", vcount(g)),
paste0("Tamaño = ", ecount(g)),
paste0("Densidad = ", round(edge_density(g), 2)),
paste0("Grado promedio = ", round(mean(degree(g)), 2)),
paste0("Diámetro = ", round(diameter(g, directed = FALSE, unconnected = TRUE), 2)),
paste0("Distancia promedio = ", round(mean_distance(g, directed = FALSE, unconnected = TRUE), 2)),
paste("Triángulos = ", sum(count_triangles(g))/2),
paste("Transitividad = ", round(transitivity(g, type="global"), 2)),
paste("Num. clique = ", clique_num(graph = g))), bty = "n", cex = 0.7, text.width = 2)
Se definen las mismas categorías de exposición de 4 niveles utilizadas previamente:
- \(c_{11}\): Exposición directa e indirecta
- \(c_{10}\): Exposición directa aislada
- \(c_{01}\): Exposición indirecta
- \(c_{00}\): Sin exposición
En la red de conflicto, se observó que 113 estudiantes no estuvieron expuestos, 125 recibieron solo exposición indirecta, 10 tuvieron exposición directa aislada, y 14 estuvieron expuestos tanto de manera directa como indirecta.
# Vector de asignación
z <- ifelse(V(g)$TREAT == "(1) Treatment", 1, 0)
# Identificadora
A <- as_adjacency_matrix(g)
I.ex.nbrs <- as.numeric(z%*%A > 0)
# Categorías de exposición
## Directa e Indirecta
V(g)[z*I.ex.nbrs==1]$exposure <- "c['11']"
## Indirecta
V(g)[(1-z)*I.ex.nbrs==1]$exposure <- "c['01']"
## Directa aislada
V(g)[z*(1-I.ex.nbrs)==1]$exposure <- "c['10']"
## No exposición
V(g)[(1-z)*(1-I.ex.nbrs)==1]$exposure <- "c['00']"
table(V(g)$exposure)##
## c['00'] c['01'] c['10'] c['11']
## 113 125 10 14A continuación, se presentan las categorías de exposición, representadas por diferentes colores. Además, se puede observar el número de eventos disciplinarios reportados para cada estudiante, que se muestra mediante el tamaño de los vértices. Se nota que los estudiantes que recibieron exposición indirecta, así como aquellos que no fueron expuestos a la intervención, tendieron a presentar un mayor número de eventos disciplinarios a lo largo del año.
# Definir colores para los valores de exposición
palette_colors <- brewer.pal(length(unique(V(g)$exposure)), "Set3")
exposure_colors <- setNames(palette_colors, unique(V(g)$exposure))
vertex_colors <- exposure_colors[V(g)$exposure]
# Visualización
plot(g, layout = l, vertex.label = NA, vertex.color = vertex_colors, vertex.size = sqrt(V(g)$event_disc)*4)
legend("topleft",
legend = parse(text = unique(V(g)$exposure)),
fill = exposure_colors,
title = "Niveles de exposición")
legend("left",
legend = as.character(sort(unique(V(g)$event_disc))),
pch = 21,
pt.bg = "grey", # Color de fondo de los puntos en la leyenda
pt.cex = sqrt(1:length(V(g)$event_disc)),
title = "Eventos \ndisciplinarios")
Las teorías del comportamiento humano sugieren que los individuos prestan atención al comportamiento de ciertas personas de su comunidad para comprender qué es socialmente normativo y ajustar su propio comportamiento en respuesta. Específicamente, el experimento realizado por Paluck et. al (2020) intenta determinar el poder relativo de ciertos individuos para influir en el comportamiento de otros, por medio de una intervención experimental para fomentar el comportamiento anticonflicto entre un grupo de estudiantes.
Teniendo en cuenta lo anterior, se espera que el comportamiento conflictivo, medido por el número de eventos disciplinarios en los que los estudiantes se involucran a lo largo del año, disminuya en función de la exposición a la intervención.
O.c00 <- mean(V(g)$event_disc[V(g)$exposure == "c['00']"])
O.c01 <- mean(V(g)$event_disc[V(g)$exposure == "c['01']"])
O.c10 <- mean(V(g)$event_disc[V(g)$exposure == "c['10']"])
O.c11 <- mean(V(g)$event_disc[V(g)$exposure == "c['11']"])
c(O.c00, O.c01, O.c10, O.c11)## [1] 0.9734513 1.0240000 0.7000000 1.5714286Sin embargo, los resultados sugieren lo contrario. En orden de mayor a menor conflictividad, en promedio, se observa lo siguiente: exposición completa, exposición indirecta, sin exposición y exposición directa aislada.
c(O.c01, O.c10, O.c11) - O.c00## [1] 0.05054867 -0.27345133 0.59797724Se pueden usar simulaciones de Monte Carlo para aproximar el valor de las probabilidades de exposición. Con solo 24 de 262 estudiantes tratados, solo las probabilidades de exposición indirecta y sin exposición tienen valores altos.
# Inicializar
set.seed(123)
nv <- vcount(g)
m <- 24 # Número de expuestos directamente
n <- 10000 # Número de ensayos de Montecarlo
# 4 condiciones de exposición
I11 <- matrix(,nrow=nv,ncol=n)
I10 <- matrix(,nrow=nv,ncol=n)
I01 <- matrix(,nrow=nv,ncol=n)
I00 <- matrix(,nrow=nv,ncol=n)
# Muestreo Monte Carlo
for(i in 1:n){
z <- rep(0, nv) # Vector z
reps.ind <- sample((1:nv), m, replace = FALSE)
z[reps.ind] <- 1 # Asignar 1 al vector z
reps.nbrs <- as.numeric(z%*%A > 0) # Indicador de vecinos expuestos
# Una matriz por cada condición de exposición
I11[,i] <- z*reps.nbrs
I10[,i] <- z*(1-reps.nbrs)
I01[,i] <- (1-z)*reps.nbrs
I00[,i] <- (1-z)*(1-reps.nbrs)
}
I11.11 <- I11%*%t(I11)/n
I10.10 <- I10%*%t(I10)/n
I01.01 <- I01%*%t(I01)/n
I00.00 <- I00%*%t(I00)/n
# Visualización
plot_Ik <- function(M, title) {
# Márgenes
par(mar = c(4,4,4,4))
# Matriz
image(M, zlim = c(0, 0.7), xaxt = "n", yaxt = "n", col = cm.colors(16))
# Título
mtext(side = 3, text = title, at = 0.5, las = 1)
# Añadir barra de colores
image.plot(legend.only = TRUE, zlim = range(seq(0.0,0.7,0.1)), col = cm.colors(16))
}
# Crear los cuatro gráficos
par(mfrow = c(2, 2))
plot_Ik(I11.11, expression("Exposición Directa + Indirecta"~(c["11"])))
plot_Ik(I10.10, expression("Exposición Directa Aislada"~(c["10"])))
plot_Ik(I01.01, expression("Exposición Indirecta"~(c["01"])))
plot_Ik(I00.00, expression("No Exposición"~(c["00"])))
Por último, las estimaciones de los efectos causales promedio \(\tau(c_{11}, c_{00})\), \(\tau(c_{10}, c_{00})\) y \(\tau(c_{01}, c_{00})\) se obtienen a continuación.
# Vector de asignación
z <- ifelse(V(g)$TREAT == "(1) Treatment", 1, 0)
# Indicador de vecinos expuestos
reps.nbrs <- as.numeric(z%*%A > 0)
# Condiciones de exposición
c11 <- z*reps.nbrs
c10 <- z*(1-reps.nbrs)
c01 <- (1-z)*reps.nbrs
c00 <- (1-z)*(1-reps.nbrs)
# Estimadores Horvitz y Thompson
Obar.c11 <- O.c11*mean(c11/diag(I11.11), na.rm = TRUE)
Obar.c10 <- O.c10*mean(c10/diag(I10.10), na.rm = TRUE)
Obar.c01 <- O.c01*mean(c01/diag(I01.01), na.rm = TRUE)
Obar.c00 <- O.c00*mean(c00/diag(I00.00), na.rm = TRUE)
# Average causal effects
ACE <- c(Obar.c11, Obar.c10, Obar.c01) - Obar.c00
ACE## [1] 1.0430773 -0.2918506 0.4255092# Diferencias con el valor obtenido
print(ACE - c(O.c11-O.c00, O.c10-O.c00, O.c01-O.c00))## [1] 0.44510005 -0.01839927 0.374960506 Referencias
Kolaczyk, E. D., & Csárdi, G. (2014). Statistical analysis of network data with R (Vol. 65). New York: Springer. Chapter 10.
Paluck, Elizabeth Levy, Shepherd, Hana R., and Aronow, Peter. Changing Climates of Conflict: A Social Network Experiment in 56 Schools, New Jersey, 2012-2013. Inter-university Consortium for Political and Social Research [distributor], 2020-09-14. https://doi.org/10.3886/ICPSR37070.v2
7 Anexo
\[ V(\hat{\tau}_{ATE}(z)) = \frac{1}{N^2}V( \sum_{i=1}^{N} \left[ \frac{z_i}{\pi_{i}}y_{i1} - \frac{(1 - z_i)}{1-\pi_{i}} y_{i0} \right] ) \] \[ = \frac{1}{N^2}V( \sum_{i=1}^{N} \left[ (\frac{y_{i1}}{\pi_{i}}+ \frac{y_{i0}}{1-\pi_{i}})z_i- \frac{y_{i0}}{1-\pi_{i}} \right] ) \] Dado que bajo un diseño experimental completamente aleatorizado la probabilidad de asignación de todas los individuos es constante, se tiene que \[ =\frac{1}{N^2}\sum_i \left[ \left( \frac{y_{i1}}{\pi} + \frac{y_{i0}}{1-\pi} \right)^2 V(Z_i) \right] + \sum_{i \ne j} \sum \left[ \left( \frac{y_{i1}}{\pi} + \frac{y_{i0}}{1-\pi} \right) \left( \frac{y_{j1}}{\pi} + \frac{y_{j0}}{1-\pi} \right) \text{Cov}(Z_i, Z_j) \right] \] \[ =\frac{1}{N^2} \sum_{i} \sum_{j} \left[ \left( \frac{y_{i1}}{\pi} + \frac{y_{i0}}{1-\pi} \right) \left( \frac{y_{j1}}{\pi} + \frac{y_{j0}}{1-\pi} \right) \text{Cov}(Z_i, Z_j) \right] \] Al denotar la covarianza \(\text{Cov}(Z_i, Z_j)\) como \(\Delta_{ij}\) y al aplicar la definición de Yates & Grundy de la varianza del estimador de Horvitz-Thompson, se obtiene \[ =\frac{1}{N^2}\left(-\frac{1}{2} \sum_{i} \sum_{j} \Delta_{ij} \left[ \left( \frac{y_{i1}}{\pi} + \frac{y_{i0}}{1-\pi} \right) - \left( \frac{y_{j1}}{\pi} + \frac{y_{j0}}{1-\pi} \right) \right]^2\right) \] Dado que \(\pi=\frac{T}{N}\) y \(1-\pi=\frac{(N-T)}{N}\), \[ =\frac{1}{N^2}\left(-\frac{1}{2} \sum_{i \ne j} \sum \Delta_{ij} \frac{N^2}{T^2\left(N-T\right)^2} \left[\left(N-T\right)\left(y_{i1} - y_{j1}\right)- T\left(y_{i0} - y_{j0}\right) \right]^2\right), (1) \] Veamos a qué es igual \(\Delta_{ij}\) para todo \(i \ne j \in V\), \[ \Delta_{ij} = \pi_{ij}-\pi_{i}\pi_{j} = \frac{T(T-1)}{N(N-1)}-\frac{T(^2}{N^2}=\frac{-T(N-T)}{N^2(N-1)} \] Sustituyendo esta expresión en (1), se obtiene \[ =\frac{1}{N^2}\left(-\frac{1}{2} \sum_{i} \sum_{j} \frac{-T(N-T)}{N^2(N-1)} \frac{N^2}{T^2\left(N-T\right)^2} \left[\left(N-T\right)\left(y_{i1} - y_{j1}\right)- T\left(y_{i0} - y_{j0}\right) \right]^2\right), \] \[ =\frac{1}{2N^2}\left(\left[\frac{1}{T(N-T)}\right]\left[\frac{1}{N-1}\right] \left[\left(N-T\right)^2\sum_{i} \sum_{j} \left(y_{i1} - y_{j1}\right)^2 - 2\left(N-T\right)T \sum_{i} \sum_{j}\left(y_{i1} - y_{j1}\right)\left(y_{i0} - y_{j0}\right) + T^2\sum_{i} \sum_{j}\left(y_{i0} - y_{j0}\right)^2\right]\right) (2) \] Desarrollando los términos de la sumatoria por separado, se tiene que \[ \sum_{i} \sum_{j} \left(y_{i1} - y_{j1}\right)^2 = \sum_{i} \sum_{j} \left( \left(y_{i1} - \overline{y_{1}}\right) - \left(y_{j1} - \overline{y_{1}}\right) \right)^2 \]
\[ =\sum_{i} \sum_{j} \left(y_{i1} - \overline{y_{1}}\right)^2 - 2\sum_{i} \sum_{j}\left(y_{i1} - \overline{y_{1}}\right) \left(y_{j1} - \overline{y_{1}}\right)+\sum_{i} \sum_{j} \left(y_{j1} - \overline{y_{1}}\right)^2 \]
Dado que \(\sum_{i}\left(y_{i1} - \overline{y_{1}}\right) \sum_{j}\left(y_{j1} - \overline{y_{1}}\right)=0\), entonces \[ =\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2\sum_{j}1 +\sum_{j} \left(y_{j1} - \overline{y_{1}}\right)^2 \sum_{i} 1 \] \[ =N\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2+N\sum_{j} \left(y_{j1} - \overline{y_{1}}\right)^2 \] \[ =2N\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 (3) \] Analógamente, se tiene que, \[ \sum_{i} \sum_{j} \left(y_{i0} - y_{j0}\right)^2 =2N\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 (4) \] \[ \sum_{i} \sum_{j}\left(y_{i1} - y_{j1}\right)\left(y_{i0} - y_{j0}\right)=2N\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)\left(y_{i0} - \overline{y_{0}}\right) (5) \] Sustituyendo (3), (4) y (5) en (2), \[ =\frac{1}{2N^2}\left[\frac{1}{T(N-T)}\right]\left[\frac{1}{N-1}\right] \left[\left(N-T\right)^2 \left(2N\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2\right) - 2\left(N-T\right)T\left(2N\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)\left(y_{i0} - \overline{y_{0}}\right)\right) + T^2\left(2N\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2\right) \right] \] \[ =\frac{1}{N}\left[\frac{1}{N-1}\right] \left[ \left[\frac{N-T}{T}\right]\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 - 2\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)\left(y_{i0} - \overline{y_{0}}\right) + \frac{T}{N-T}\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 \right] \]
\[ =\frac{1}{N}\left[\frac{1}{N-1}\right] \left[ \left[\frac{N-T}{T}\right]\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 - 2\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)\left(y_{i0} - \overline{y_{0}}\right) + \frac{T}{N-T}\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 \right] \]
Sumando y restando , se tiene que \[ =\frac{1}{N}\left[\frac{1}{N-1}\right] \left[ \left[\frac{C-T}{T}\right]\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 + \sum_{i}\left((y_{i1} - \overline{y_{1}}) -(y_{i0} - \overline{y_{0}})\right)^2 + \frac{T-C}{C}\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 \right] \] \[ =\frac{1}{N}\left[\frac{1}{N-1}\right] \left[ \left[\frac{C-T}{T}\right]\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 + \sum_{i}\left((y_{i1} - y_{i0}) - (\overline{y_{1}} - \overline{y_{0}})\right)^2 + \frac{T-C}{C}\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 \right] \] \[ =\frac{1}{N}\left[\frac{1}{N-1}\right] \left[\frac{C-T}{T}\sum_{i}\left(y_{i1} - \overline{y_{1}}\right)^2 + \frac{T-C}{C}\sum_{i}\left(y_{i0} - \overline{y_{0}}\right)^2 + \sum_{i}\left((y_{i1} - y_{i0}) - \tau_{ATE}\right)^2 \right] \]